% Natasha Savage



function [rt,m,em,ec,rd,rdim,rdic,ic,RDic,RDIic,recic,recm,reca] = Euler_step_GDI42ic(RTn,Mn,Emn,Ecn,RDn,RDImn,RDIcn,In,RDicn,RDIicn,RecCmplxicn,RecCmplxn,RecAn,h)


global k1a 
global k1b 
global k2a  
global k2b 
global k3 
global k4a 
global k4b 
global k5a 
global k5b 
global k6a 
global k6b 
global k7 
global eta
global dx
global dx2
global kbirth %Masha dev sept2012
global k8a %rate of pheromone binding to receptor complex                                   %Masha dev sept2012
global k8b %rate of pheromone unbinding to receptor complex                                 %Masha dev sept2012
global k9 %rate of nucleotide exchange in Cdc42D->T assisted by liganded receptor RecA %Masha dev sept2012

rt    = RTn    + h.* ((k2a*Emn+k3*Mn+k9*RecAn).*RDn - k2b*RTn - k4a*Emn.*RTn + k4b*Mn - k7*Ecn.*RTn);   

m     = Mn     + h.*(k4a*Emn.*RTn - k4b*Mn + k7*Ecn.*RTn);             

em    = Emn    + h.*(k1a*Ecn - k1b*Emn - k4a*Emn.*RTn + k4b*Mn);

ec    = Ecn    + h.*(eta*(k1b*Emn - (k1a+k7*RTn).*Ecn));

rd    = RDn    + h.*(k2b*RTn - (k2a*Emn+k3*Mn+k9*RecAn).*RDn + k6b*RDImn -k6a*In.*RDn); % Added GDP->GTP  exchange by Reca (k9) %Masha dev sept2012

rdim  = RDImn  + h.*(k6a*In.*RDn - k6b*RDImn + k5a*RDIcn - k5b*RDImn);

rdic  = RDIcn  + h.*(eta*(k5b*RDImn  - k5a*RDIcn) + ...
                     eta*(k5b*RDIicn - k5a*mean(mean(RDIcn)))*dx2*dx2/dx/dx );

ic    = In     + h.*(eta*(k6b*RDImn  - k6a*In.*RDn) + ...
                     eta*(k6b*RDIicn - k6a*mean(mean(In))*RDicn)*dx2*dx2/dx/dx ); 

% internal compartment concentration is a scalar
% so use the mean of the cytoplasm
RDic  = RDicn  + h.*(  k6b*RDIicn - k6a*mean(mean(In))*RDicn);
RDIic = RDIicn + h.*(- k6b*RDIicn + k6a*mean(mean(In))*RDicn + k5a*mean(mean(RDIcn)) - k5b*RDIicn);

%Masha dev sept2012
recic = RecCmplxicn + h.*(kbirth); %receptor secretion in the internal compartment
recm  = RecCmplxn   + h.*(-k8a*RecCmplxn + k8b*RecAn); %free receptor on the membrane convertion to being bound to pheromone
reca  = RecAn       + h.*( k8a*RecCmplxn - k8b*RecAn); %receptor bound to pheromone time evolution

end